GBIF_ENM_Ensambles
Velasco-Vinasco J.A. y Rodriguez-Rodriguez J.D.
2024-10-28
Descarga de datos de GBIF
Este ejercicio estará enfocado en descargar datos del repositorio GBIF y usar algunas herramientas sencillas para la limpieza de datos.
##Instalar los paquetes
#install.packages("sf")
#install.packages("sp")
#install.packages("raster")
#install.packages("dismo")
#install.packages("spThin")
#install.packages("ecospat")
#install.packages("geodata")
#install.packages("sdm")
#install.packages("usdm")
#install.packages("biomod2")
#install.packages("rasterVis")
#Activar los paquetes
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUElibrary(sp)
library(raster)
library(dismo)
library(spThin)## Loading required package: spam## Spam version 2.10-0 (2023-10-23) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: grid## Loading required package: fields## Loading required package: viridisLite##
## Try help(fields) to get started.## Loading required package: knitrlibrary(ecospat)
library(geodata)## Loading required package: terra## terra 1.7.83##
## Attaching package: 'terra'## The following object is masked from 'package:knitr':
##
## spin## The following object is masked from 'package:fields':
##
## describe## The following object is masked from 'package:grid':
##
## depth##
## Attaching package: 'geodata'## The following object is masked from 'package:fields':
##
## worldlibrary(sdm)## sdm 1.2-46 (2024-07-17)library(usdm)
library(rgdal)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:2.1-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(biomod2)## biomod2 4.2-5-2 loaded.
## /!\ New set up for modeling options. We apologize for the trouble ^[*.*]^## Loading required package: nnet## Loading required package: rpart## Loading required package: mda## Loading required package: class## Loaded mda 0.5-4## Loading required package: gam## Loading required package: splines## Loading required package: foreach## Loaded gam 1.22-5## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:usdm':
##
## Variogram## The following object is masked from 'package:raster':
##
## getData## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.##
## Attaching package: 'mgcv'## The following objects are masked from 'package:gam':
##
## gam, gam.control, gam.fit, s## The following object is masked from 'package:nnet':
##
## multinom## Loading required package: gbm## Loaded gbm 2.2.2## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3## Loading required package: earth## Loading required package: Formula## Loading required package: plotmo## Loading required package: plotrix##
## Attaching package: 'plotrix'## The following object is masked from 'package:terra':
##
## rescale## The following object is masked from 'package:fields':
##
## color.scale## Loading required package: maxnet## Loading required package: randomForest## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.## Loading required package: xgboostlibrary(rasterVis)## Loading required package: lattice##
## Attaching package: 'rasterVis'## The following object is masked from 'package:gam':
##
## gplot#Establacemos el directorio de trabajo
setwd("/home2/taller_R/")
getwd()## [1] "/home2/taller_R"Buscaremos datos de presencia para una especie de lagartija del desierto (Phrynosoma modestum)
sp.gbif <- dismo::gbif(genus="Phrynosoma", species="modestum", geo=TRUE, removeZeros=TRUE, download=TRUE)## 5012 records found## 0-300-600-900-1200-1500-1800-2100-2400-2700-3000-3300-3600-3900-4200-4500-4800-5012 records downloadedVamos a indexar solos los datos de las coordenadas geográficas. Debemos saber en que posición están.
names(sp.gbif)## [1] "acceptedNameUsage" "acceptedScientificName"
## [3] "acceptedTaxonKey" "accessRights"
## [5] "adm1" "adm2"
## [7] "associatedOccurrences" "associatedReferences"
## [9] "associatedSequences" "associatedTaxa"
## [11] "basisOfRecord" "bibliographicCitation"
## [13] "catalogNumber" "class"
## [15] "classKey" "cloc"
## [17] "collectionCode" "collectionID"
## [19] "collectionKey" "continent"
## [21] "coordinatePrecision" "coordinateUncertaintyInMeters"
## [23] "country" "crawlId"
## [25] "datasetID" "datasetKey"
## [27] "datasetName" "dateIdentified"
## [29] "day" "depth"
## [31] "depthAccuracy" "disposition"
## [33] "distanceFromCentroidInMeters" "dynamicProperties"
## [35] "earliestAgeOrLowestStage" "earliestEonOrLowestEonothem"
## [37] "earliestEpochOrLowestSeries" "earliestPeriodOrLowestSystem"
## [39] "elevation" "elevationAccuracy"
## [41] "endDayOfYear" "establishmentMeans"
## [43] "eventDate" "eventRemarks"
## [45] "eventTime" "family"
## [47] "familyKey" "fieldNotes"
## [49] "fieldNumber" "footprintSRS"
## [51] "footprintWKT" "formation"
## [53] "fullCountry" "gbifID"
## [55] "gbifRegion" "genericName"
## [57] "genus" "genusKey"
## [59] "geodeticDatum" "georeferencedBy"
## [61] "georeferencedDate" "georeferenceProtocol"
## [63] "georeferenceRemarks" "georeferenceSources"
## [65] "georeferenceVerificationStatus" "habitat"
## [67] "higherClassification" "higherGeography"
## [69] "higherGeographyID" "hostingOrganizationKey"
## [71] "http://unknown.org/captive" "http://unknown.org/nick"
## [73] "http://unknown.org/recordID" "identificationID"
## [75] "identificationQualifier" "identificationRemarks"
## [77] "identifiedBy" "identifier"
## [79] "individualCount" "informationWithheld"
## [81] "installationKey" "institutionCode"
## [83] "institutionID" "institutionKey"
## [85] "isInCluster" "ISO2"
## [87] "isSequenced" "iucnRedListCategory"
## [89] "key" "kingdom"
## [91] "kingdomKey" "language"
## [93] "lastCrawled" "lastInterpreted"
## [95] "lastParsed" "lat"
## [97] "latestAgeOrHighestStage" "latestEpochOrHighestSeries"
## [99] "latestPeriodOrHighestSystem" "license"
## [101] "lifeStage" "locality"
## [103] "locationAccordingTo" "locationID"
## [105] "locationRemarks" "lon"
## [107] "materialSampleID" "modified"
## [109] "month" "municipality"
## [111] "nameAccordingTo" "namePublishedIn"
## [113] "namePublishedInYear" "nomenclaturalCode"
## [115] "occurrenceID" "occurrenceRemarks"
## [117] "occurrenceStatus" "organismID"
## [119] "organismName" "organismQuantity"
## [121] "organismQuantityType" "otherCatalogNumbers"
## [123] "ownerInstitutionCode" "parentNameUsage"
## [125] "phylum" "phylumKey"
## [127] "preparations" "previousIdentifications"
## [129] "projectId" "protocol"
## [131] "publishedByGbifRegion" "publishingCountry"
## [133] "publishingOrgKey" "recordedBy"
## [135] "recordNumber" "references"
## [137] "rights" "rightsHolder"
## [139] "samplingProtocol" "scientificName"
## [141] "scientificNameID" "sex"
## [143] "species" "speciesKey"
## [145] "specificEpithet" "startDayOfYear"
## [147] "taxonConceptID" "taxonID"
## [149] "taxonKey" "taxonomicStatus"
## [151] "taxonRank" "taxonRemarks"
## [153] "type" "typeStatus"
## [155] "typifiedName" "verbatimCoordinateSystem"
## [157] "verbatimElevation" "verbatimEventDate"
## [159] "verbatimLocality" "verbatimSRS"
## [161] "vernacularName" "year"
## [163] "downloadDate"sp.coords <- na.omit(sp.gbif[,c(106, 96, 140)]) # lon, lat, species
head(sp.coords)## lon lat scientificName
## 1 -105.4433 32.23785 Phrynosoma modestum Girard, 1852
## 2 -105.4444 32.23727 Phrynosoma modestum Girard, 1852
## 3 -106.5001 32.86959 Phrynosoma modestum Girard, 1852
## 4 -107.4712 32.96034 Phrynosoma modestum Girard, 1852
## 5 -103.6514 29.57966 Phrynosoma modestum Girard, 1852
## 6 -106.7387 32.38824 Phrynosoma modestum Girard, 1852Graficamos los puntos de las observaciones
plot(sp.coords$lon, sp.coords$lat) 
Aquí vemos que hay muchos puntos muy agrupados. Este tipo de datos puede
influir negativamente en la calibración del modelo de nicho y la
predicción geográfica. Vamos a tratar de reducir ese efecto adelgazando
la muestra de puntos. Lo primero que hacemos es tratar de eliminar
puntos a una distancia menor de 25 km. Como tenemos un tamaño de muestra
grande (> 3000 registros) podemos hacer esto sin mucho problema.
Usaremos el paquete spThin para esto.
sp.thin <-thin(sp.coords, lat.col="lat", long.col="lon", spec.col="scientificName",
thin.par=25, reps=1, locs.thinned.list.return=TRUE, write.files=TRUE,
max.files=1, out.dir=".", out.base = "modestum.train.25km", write.log.file=TRUE, log.file="modestum.log.txt",
verbose=TRUE)## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Mon Oct 28 15:18:34 2024## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : There appear to be more than one species name in this *.csv
## file.## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : Only using species name: Phrynosoma modestum Girard, 1852##
## lat.long.thin.count
## 455
## 1
## [1] "Maximum number of records after thinning: 455"
## [1] "Number of data.frames with max records: 1"
## [1] "Writing new *.csv files"## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : Writing new *.csv files to output directory: .## [1] "Writing file: ./modestum.train.25km_thin1.csv"Ahora vamos primero a identificar que tipo de objeto se acaba de generar:
class(sp.thin)## [1] "list"str(sp.thin)## List of 1
## $ :'data.frame': 455 obs. of 2 variables:
## ..$ Longitude: num [1:455] -107 -108 -104 -105 -109 ...
## ..$ Latitude : num [1:455] 32.9 32.6 30.2 35.3 32.3 ...head(sp.thin[[1]])## Longitude Latitude
## 3 -106.5001 32.86959
## 8 -107.9779 32.56886
## 26 -103.5483 30.15385
## 31 -105.2019 35.30268
## 33 -109.1064 32.31848
## 56 -108.8535 31.83653Tenemos que tener en cuenta como vienen organizados los datos para luego manipularlos como queramos. Vamos a graficar los puntos adelgazados a una distancia no menor de 25 km.
plot(sp.thin[[1]]$Longitude, sp.thin[[1]]$Latitude)
Si comparamos con el objeto que contenia todos los puntos podemos ver
que hemos reducido el sobre-muestreo en algunas regiones y esto nos
facilitará la calibración y posterior predicción de nuestro modelo de
nicho. Vamos a ver donde caen nuestros puntos dentro de un mapa.
Vamos a convertir nuestros puntos en un objeto espacial.
coords <- sp.thin[[1]]
coordinates(coords)<- ~Longitude+Latitude
plot(coords)
Si solo nos interesan los puntos que se distribuyen en México, podríamos
excluir los puntos que caen por fuera del país. Vamos a hacerlo. Primero
necesitamos cargar el mapa de México.
mx <- readOGR(dsn="./input/Poligono_Mexico", layer="gadm36_MEX_1")## Warning in readOGR(dsn = "./input/Poligono_Mexico", layer = "gadm36_MEX_1"):
## OGR support is provided by the sf and terra packages among others## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "/home2/taller_R/input/Poligono_Mexico", layer: "gadm36_MEX_1"
## with 32 features
## It has 10 fieldsplot(mx)
Ahora le ponemos los puntos encima al mapa
plot(coords)
plot(mx, add=TRUE)
Con un clip vamos a quitar los puntos que caen por fuera. Pero primero
definimos el sistema de coordenadas para nuestros puntos
mx## class : SpatialPolygonsDataFrame
## features : 32
## extent : -118.3689, -86.71014, 14.53292, 32.71804 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : GID_0, NAME_0, GID_1, NAME_1, VARNAME_1, NL_NAME_1, TYPE_1, ENGTYPE_1, CC_1, HASC_1
## min values : MEX, Mexico, MEX.1_1, Aguascalientes, NA, NA, Distrito Federal, Federal District, NA, MX.AG
## max values : MEX, Mexico, MEX.9_1, Zacatecas, NA, NA, Estado, State, NA, MX.ZAcrdref <- CRS("+proj=longlat +datum=WGS84 +no_defs")
projection(coords)<-crdref
coords## class : SpatialPoints
## features : 455
## extent : -118.1731, -95.86836, 21.20054, 44.71625 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defscoords.clip <- coords[mx, ]
coords.clip## class : SpatialPoints
## features : 195
## extent : -112.6083, -97.70008, 21.20054, 31.775 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defsGraficamos
plot(coords.clip)
plot(mx, add=TRUE)
Ahora tenemos nuestros puntos únicamente para México y los vamos a
convertir en un objeto tipo data.frame con el fin de guardarlos en
formato CSV y explorarlos con excel.
coords.clipDF <- as.data.frame(coords.clip)
head(coords.clipDF)## Longitude Latitude
## 119 -101.6513 23.88682
## 143 -101.2116 29.33872
## 187 -101.0510 22.76873
## 201 -102.7829 28.57397
## 309 -107.1414 31.54090
## 324 -102.0906 23.06960write.csv(coords.clipDF, file="./output/Coordenadas_GBIF_Mexico.csv")
coordsMx <- read.csv("./output/Coordenadas_GBIF_Mexico.csv")
head(coordsMx)## X Longitude Latitude
## 1 119 -101.6513 23.88682
## 2 143 -101.2116 29.33872
## 3 187 -101.0510 22.76873
## 4 201 -102.7829 28.57397
## 5 309 -107.1414 31.54090
## 6 324 -102.0906 23.06960Modelos de nicho ecológico o de distribución geográfica
En este parte del ejercicio vamos a realizar un modelo de nicho ecológico sencillo usando el algoritmo bioclim.
coordsMx <- read.csv("./output/Coordenadas_GBIF_Mexico.csv")
head(coordsMx)## X Longitude Latitude
## 1 119 -101.6513 23.88682
## 2 143 -101.2116 29.33872
## 3 187 -101.0510 22.76873
## 4 201 -102.7829 28.57397
## 5 309 -107.1414 31.54090
## 6 324 -102.0906 23.06960Le quitamos la primera columna, no la necesitamos
coordsMx <- coordsMx[,c(2:3)]
head(coordsMx)## Longitude Latitude
## 1 -101.6513 23.88682
## 2 -101.2116 29.33872
## 3 -101.0510 22.76873
## 4 -102.7829 28.57397
## 5 -107.1414 31.54090
## 6 -102.0906 23.06960Subimos las capas climáticas
clim <- worldclim_global(var="bio",res=2.5,path="./")Gráficamos la primera
plot(clim[[1]])
clim<-stack(clim)Vamos a hacer un buffer a nuestros puntos. Este buffer será nuestra área M del Diagrama BAM y es el área que la especie hipóteticamente ha podido “explorar” durante su historia evolutiva (i.e., dispersión, colonización, extinción):
coordinates(coordsMx) <- ~Longitude+Latitude
coordsMx <- sf::st_as_sf(coordsMx)
buffer.sp <- st_buffer(coordsMx, byid=FALSE, id=NULL, dist =3.0, nQuadsegs=5,
endCapStyle="ROUND", joinStyle="ROUND", mitreLimit=1.0,dissolve = TRUE)
buffer.sp <- st_union(buffer.sp)Graficamos
plot(buffer.sp)
plot(coordsMx, add=TRUE)
plot(mx,add=TRUE)
Vamos a cortar las capas climáticas usando este buffer
buffer.sp<-as_Spatial(buffer.sp)
class(buffer.sp)## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"clim.m <- crop(clim, buffer.sp)
plot(clim.m[[1]])
Ahora vamos a calibrar el modelo usando el algoritmo Bioclim
# primero convertimos las coordenadas a dataframe
coordsDF <- st_coordinates(coordsMx)
coordsDF <- as.data.frame(coordsDF)
class(coordsDF)## [1] "data.frame"# ajustamos un modelo de bioclim a los datos con las capas climáticas
bioc.modesta <- bioclim(clim.m, coordsDF)Vamos a ver las gráficas de las variables de respuesta
response(bioc.modesta, var="wc2.1_2.5m_bio_1")
response(bioc.modesta, var="wc2.1_2.5m_bio_12")
response(bioc.modesta, var="wc2.1_2.5m_bio_4")
plot(bioc.modesta, a=1, b=2, p=0.95)
Ahora vamos a generar la predicción geográfica:
pred.bioc <- predict(clim.m, bioc.modesta)
plot(pred.bioc)
Vamos a generar una predicción binaria (0 ó 1; presencia vs. ausencia).
Usaremos el criterio de umbral del “minimum presence training”
coords.bioc <- extract(pred.bioc, coords)
head(coords.bioc)## [1] 0.02051282 0.02051282 0.01025641 NA 0.03589744 0.03076923a <- min(coords.bioc)
a## [1] NAUsaremos ese valor para generar la predicción binaria. Por encima de 0.0073 será presencia y por debajo, ausencia.
bin.bioc <- reclassify(pred.bioc,c(-Inf,a,0, a,Inf,1))
bin.bioc## class : RasterLayer
## dimensions : 398, 502, 199796 (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -115.625, -94.70833, 18.20833, 34.79167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 0.6410256 (min, max)Graficamos
plot(bin.bioc)
Vamos a usar los registros completos de la distribución de la especie para ver si tenemos un problema de truncación de nicho
Hacemos otro buffer pero con todos los puntos adelgazados. Checamos la anchura para que todo quede conectado
coords<-as.data.frame(coords)
head(coords)## Longitude Latitude
## 3 -106.5001 32.86959
## 8 -107.9779 32.56886
## 26 -103.5483 30.15385
## 31 -105.2019 35.30268
## 33 -109.1064 32.31848
## 56 -108.8535 31.83653write.csv(coords, file="./output/Coordenadas_GBIF.csv")
coords <- read.csv("./output/Coordenadas_GBIF.csv")
head(coords)## X Longitude Latitude
## 1 3 -106.5001 32.86959
## 2 8 -107.9779 32.56886
## 3 26 -103.5483 30.15385
## 4 31 -105.2019 35.30268
## 5 33 -109.1064 32.31848
## 6 56 -108.8535 31.83653coords <- coords[,c(2:3)]
head(coords)## Longitude Latitude
## 1 -106.5001 32.86959
## 2 -107.9779 32.56886
## 3 -103.5483 30.15385
## 4 -105.2019 35.30268
## 5 -109.1064 32.31848
## 6 -108.8535 31.83653coordinates(coords) <- ~Longitude+Latitude
head(coords)## class : SpatialPoints
## features : 1
## extent : -106.5001, -106.5001, 32.86959, 32.86959 (xmin, xmax, ymin, ymax)
## crs : NAcoords <- sf::st_as_sf(coords)
buffer.sp2 <- st_buffer(coords, byid=FALSE, id=NULL, dist =6.0, nQuadsegs=5,
endCapStyle="ROUND", joinStyle="ROUND", mitreLimit=1.0,dissolve = TRUE)
buffer.sp2 <- st_union(buffer.sp2)Graficamos
plot(buffer.sp2)
plot(coords, add=TRUE)
plot(mx, add=TRUE)
buffer.sp2<-as_Spatial(buffer.sp2)
class(buffer.sp2)## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"clim.m2 <- crop(clim, buffer.sp2)
plot(clim.m2[[1]])
coordsDF2 <- st_coordinates(coords)
coordsDF2 <- as.data.frame(coordsDF2)
class(coordsDF2)## [1] "data.frame"bioc.modesta.v2 <- bioclim(clim.m2, coordsDF2)
pred.bioc.v2 <- predict(clim.m2, bioc.modesta.v2)
plot(pred.bioc.v2)
Comparemos las dos predicciones:
par(mfrow=c(1,2))
plot(pred.bioc, main="M puntos en México")
plot(pred.bioc.v2, main="M todos los puntos")
Ajustar las extensiones para ambos rasters
pred.bioc <- projectRaster(from=pred.bioc, to=pred.bioc.v2, crs=crs(pred.bioc.v2), method="bilinear")
par(mfrow=c(1,2))
plot(pred.bioc, main="M puntos en México")
plot(pred.bioc.v2, main="M todos los puntos")
Aquí hemos calibrado un modelo de nicho ecológico usando todos los datos disponibles y por lo tanto no tenemos forma de evaluar que tanto error puede tener nuestra predicción geográfica. Vamos a partir los datos en dos conjuntos. Uno de calibración y uno de validación:
La partición de los datos será aleatoria. Usaremos la función kfold para esto. Usaremos 5 grupos
group <- kfold(coordsDF2, 5)
group## [1] 3 1 1 3 4 2 5 1 5 2 1 3 5 4 5 2 4 2 1 4 1 4 4 2 2 2 3 4 4 1 4 2 1 3 5 5 3
## [38] 2 3 1 3 4 3 3 1 3 4 4 4 3 1 2 4 3 2 4 2 5 4 3 3 1 3 4 4 5 1 4 1 3 5 5 5 5
## [75] 3 3 3 5 5 5 2 2 2 4 4 5 1 1 3 4 2 3 3 1 2 3 5 2 5 3 5 5 2 5 2 1 1 5 4 3 3
## [112] 5 5 3 1 1 2 5 1 4 4 2 2 5 2 1 1 1 1 5 3 4 2 2 2 4 4 4 2 2 4 3 3 1 3 4 2 2
## [149] 2 2 3 2 1 3 2 5 1 4 1 2 1 4 1 3 2 4 5 4 5 3 4 2 4 5 3 1 3 5 3 5 1 5 2 5 4
## [186] 5 4 2 4 4 3 4 5 4 5 3 1 3 1 3 4 1 1 3 2 5 5 3 3 4 4 1 5 2 2 4 2 3 5 3 3 3
## [223] 1 1 3 5 5 4 2 5 2 1 3 2 4 5 2 3 1 1 1 4 2 2 3 2 4 3 2 1 4 3 4 3 2 2 3 4 3
## [260] 5 1 3 1 1 1 4 2 2 1 1 4 4 4 4 5 4 1 4 1 5 4 1 1 4 4 4 3 5 2 5 2 4 2 3 5 2
## [297] 2 3 5 1 4 3 5 1 4 2 5 2 3 4 2 5 1 5 4 2 3 2 3 1 1 5 5 5 3 4 2 1 5 1 5 1 2
## [334] 2 3 2 2 5 5 1 5 1 2 5 3 3 3 2 5 2 5 1 1 1 5 1 2 1 1 3 4 5 2 3 1 2 3 5 5 2
## [371] 2 3 4 5 1 5 4 2 3 1 4 4 2 1 1 4 5 4 3 5 3 5 1 5 3 4 3 4 5 3 3 1 2 5 4 2 5
## [408] 1 1 4 3 4 3 4 1 5 5 5 1 2 3 2 4 2 3 3 2 2 5 4 1 1 2 1 5 1 4 4 1 2 3 5 2 5
## [445] 4 3 3 5 4 1 5 1 1 5 4length(group)## [1] 455pres_train <- coordsDF[group != 1, ]
pres_test <- coordsDF[group == 1, ]
head(pres_train)## X Y
## 1 -101.6513 23.88682
## 4 -102.7829 28.57397
## 5 -107.1414 31.54090
## 6 -102.0906 23.06960
## 7 -101.5608 24.71512
## 9 -101.1581 25.28576head(pres_test)## X Y
## 2 -101.2116 29.33872
## 3 -101.0510 22.76873
## 8 -100.8064 22.63093
## 11 -106.2099 27.87126
## 19 -100.5416 26.22126
## 21 -100.7210 25.91715Ya tenemos las presencias particionadas en dos conjuntos (calibración o “training” y validación o “testing”). Para validar nuestros modelos necesitamos las “ausencias”. No tenemos ausencias, así que vamos a generar unas ausencias aleatorias a partir de la geografía: a este tipo de ausencias se les conoce como “pseudo-ausencias”
ext <- extent(clim.m2)
backg <- randomPoints(clim.m2, n=1000, ext=ext)
colnames(backg) = c('Longitude', 'Latitude')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]Vamos a graficar los puntos
plot(clim.m2[[1]])
points(backg_train, pch='-', cex=1, col='yellow')
points(backg_test, pch='-', cex=1, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')
Ahora vamos a calibrar con el primer conjunto que extrajimos,
pres_train:
bc <- bioclim(clim.m2, pres_train)
plot(bc, a=1, b=2, p=0.85)
Vamos a validar nuestro modelo de Bioclim
# Usaremos la función evaluate de dismo.
eval.modesta <- evaluate(pres_test, backg_test, bc, clim.m2)
eval.modesta## class : ModelEvaluation
## n presences : 34
## n absences : 200
## AUC : 0.8564706
## cor : 0.5310529
## max TPR+TNR at : 0.01853354Vamos a calcular un umbral
tr <- threshold(eval.modesta, 'spec_sens')
tr## [1] 0.01853354Vamos a hacer la predicción en la geografía:
pred.bioc.v3 <- predict(clim.m2, bc, ext=ext, progress='')
pred.bioc.v3## class : RasterLayer
## dimensions : 852, 823, 701196 (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -124.1667, -89.875, 15.20833, 50.70833 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 0.621118 (min, max)Graficamos los dos mapas de predicciones, el continuo y el discreto usando el valor de umbral que calculamos antes.
par(mfrow=c(1,2))
plot(pred.bioc.v3, main='Bioclim, valores crudos')
plot(pred.bioc.v3 > tr, main='presencia/ausencia')
points(pres_train, pch='+')
Ahora vamos a generar modelos con varios algoritmos y crear un ensamble. Usaremos un conjunto de datos de un ratón (Neotoma cinerea “bushy-tailed woodrat”) de la familia Cricetidae que se distribuye en Estados Unidos y Cánada. Usaremos el paquete de R “sdm” (Naimi & Aráujo 2016; https://github.com/babaknaimi/sdm) que es bastante versátil.
datos <- read.csv("Neotoma.csv")
plot(datos$decimalLongitude,datos$decimalLatitude)
bios <- stack(list.files(path="./ClimaActual/", pattern = "*.tif", full.names = T))
plot(bios[[1]])
pts <- datos[,c(2,3)]
selected <- sample(1:nrow(pts), nrow(pts) * 0.5)
train <- pts[selected,] # this is the selection to be used for model training
test <- pts[-selected,] # this is the opposite of the selection which will be used for model testing
pts.all <- rbind(train,test)
# convertir el primer raster de los bios en puntos
bios.pts<-rasterToPoints(bios[[1]],fun=NULL, spatial=TRUE)
bios.all<-raster::extract(bios[[1]],bios.pts)
bios.all<-data.frame(coordinates(bios.pts),bios.all)
# generar pseudo-ausencias para los datos de calibración (50,000 o más)
presences_train<- train
dim(presences_train)## [1] 216 2pseudo_abs_train<- ecospat.rand.pseudoabsences(nbabsences=1000,
glob=bios.all,colxyglob=1:2, colvar=1:2, presence=presences_train,
colxypresence=1:2, mindist=0.1)
pseudo_abs_train<- pseudo_abs_train[,1:2]
presences_train["Rep1"] <- 1
presences_train <- plyr::rename(presences_train, replace=c("decimalLongitude"="x", "decimalLatitude"="y"))
pseudo_abs_train["Rep1"] <- 0
# Combino las presencias y las pseudo-absences del conjunto de calibración
DataSpeciesTrain.spp <-rbind(presences_train, pseudo_abs_train)
# generar pseudo-ausencias para los datos de validación
presences_test <- test
presences_test["Rep1"] <- 1
presences_test <- plyr::rename(presences_test, replace=c("decimalLongitude"="x", "decimalLatitude"="y"))
dim(presences_test)## [1] 217 3pseudo_abs_test<- ecospat.rand.pseudoabsences(nbabsences=1000,
glob=bios.all,colxyglob=1:2, colvar=1:2, presence=presences_test,
colxypresence=1:2, mindist=0.1)
pseudo_abs_test<-pseudo_abs_test[,1:2]
pseudo_abs_test["Rep1"] <- 0
# Combino las presencias y las pseudo-absences del conjunto de validación
DataSpeciesTest.spp <- rbind(presences_test, pseudo_abs_test)
# Convierto los dataframes en objectos espaciales
DataSpecies <- SpatialPointsDataFrame(coords =DataSpeciesTrain.spp[,c(1,2)], data = DataSpeciesTrain.spp,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
TestSpecies <- SpatialPointsDataFrame(coords =DataSpeciesTest.spp[,c(1,2)], data = DataSpeciesTest.spp,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
DataSpecies <- DataSpecies[,c(-1,-2)]
TestSpecies <- TestSpecies[,c(-1,-2)]Voy a examinar la colinealidad entre las variables usando el VIF (https://en.wikipedia.org/wiki/Variance_inflation_factor)
biosDF <- as.data.frame(bios)
# La función vifcor identifica las variables que deben ser excluidas
v1 <- usdm::vifcor(biosDF, th=0.7)
bios.red <- usdm::exclude(bios, v1)
bios.red## class : RasterStack
## dimensions : 1363, 2741, 3735983, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -166.8333, -52.625, 14.54167, 71.33333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : bio_05, bio_12, bio_15
## min values : -27, 35, 5
## max values : 456, 5152, 138Vamos a ajustar los datos a los requirimientos específicos del paquete sdm usando todas las bios
d1 <- sdmData(formula=Rep1~., train=DataSpecies, test=TestSpecies, predictors=bios)
d1## class : sdmdata
## ===========================================================
## number of species : 1
## species names : Rep1
## number of features : 5
## feature names : bio_01, bio_05, bio_06, ...
## type : Presence-Absence
## has independet test data? : TRUE
## number of records : train-> 1214; test-> 1216
## has Coordinates? : TRUEEstrategia de submuestreo: 70%-30% (calibración-validación) con 2 (haga 10 o más en su casa) replicas de bootstrapping. Vamos a probar maxent con default settings (ya vimos en clase qué consecuencias trae eso), random forest (los mejores en mi opinión), glm y bioclim. Se pueden ajustar más algoritmos. Vaya a la ayuda de la función tecleando ?sdm para ver más opciones.
m1 <- sdm(Rep1~., data=d1, methods=c('rf','glm','bioclim'),
replicatin='sub', test.percent=30, n=2)
m1## class : sdmModels
## ========================================================
## number of species : 1
## number of modelling methods : 3
## names of modelling methods : rf, glm, bioclim
## replicate.methods (data partitioning) : subsampling
## number of replicates (each method) : 2
## toral number of replicates per model : 2 (per species)
## test percentage (in subsampling) : 30
## ------------------------------------------
## model run success percentage (per species) :
## ------------------------------------------
## method Rep1
## ----------------------
## rf : 100 %
## glm : 100 %
## bioclim : 100 %
##
## ###################################################################
## model Mean performance (per species), using independent test dataset:
## -------------------------------------------------------------------------------
##
## ## species : Rep1
## =========================
##
## methods : AUC | COR | TSS | Deviance
## -------------------------------------------------------------------------
## rf : 0.92 | 0.65 | 0.75 | 0.53
## glm : 0.88 | 0.56 | 0.62 | 0.64
## bioclim : 0.81 | 0.52 | 0.6 | 0.96Vamos a ver las curvas ROC que evaluan el balance entre la especificidad y sensibilidad. Recuerden la clase de la curva ROC y matriz de confusión.
roc(m1)
Generar un ensamble con los mejores modelos usando el método de promedio
ponderado (weighted averaging) con base en la métrica (TSS: True
Statistics Skill)
e1 <- ensemble(m1, newdata=bios, filename='e1.img', setting=list(method='weighted',stat='TSS'))## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.Graficamos el modelo
plot(e1)## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.## Warning in x@ptr$readStart(): GDAL Message 1: NaN converted to INT_MAX.
Generamos un mapa binario usando un criterio de umbral de mínima área
predicha y donde le asignamos un error de 10% a los datos de presencia
que no usamos para calibrar.
xy <- train
coordinates(xy) <- ~decimalLongitude+decimalLatitude
mpa.e1 <- ecospat.mpa(e1, xy, perc = .9)## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.## Warning in x@ptr$extractCell(i - 1): GDAL Message 1: NaN converted to INT_MAX.mpa.e1## 10%
## 0.301##Cambiar e1 de rasterlayer a SpatRaster
e1<-rast(e1)## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.bin.e1 <- ecospat.binary.model(e1, mpa.e1)## Warning in x@ptr$classify(as.vector(rcl), NCOL(rcl), right, include.lowest, :
## GDAL Message 1: NaN converted to INT_MAX.plot(bin.e1)
Vamos a generar unas proyecciones a futuro.
# Subir bioclimáticos de escenarios futuros (RCP8.5 2070)
fut.bios <- stack(list.files(path="./CMIP5", pattern = "tif", full.names = T))
fut.bios <- projectRaster(from=fut.bios, to=bios[[1]], crs=fut.bios, method="bilinear")
# recortar las bios a la extensión
fut.bios.m <- raster::crop(fut.bios, bios[[1]])
# generar la proyección futura sobre el ensamble
e1.fut <- ensemble(m1,newdata=fut.bios.m,filename='e1_fut.img',overwrite = T, setting=list(method='weighted',stat='TSS'))## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.Graficamos
plot(e1.fut)## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.## Warning in x@ptr$readStart(): GDAL Message 1: NaN converted to INT_MAX.
Generar el binario para el escenario futuro
e1.fut<-rast(e1.fut)## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.bin.e1.fut <- ecospat.binary.model(e1.fut, mpa.e1)## Warning in x@ptr$classify(as.vector(rcl), NCOL(rcl), right, include.lowest, :
## GDAL Message 1: NaN converted to INT_MAX.Graficamos el mapa binario proyectado a futuro
plot(bin.e1.fut)
Ahora vamos a identificar áreas de estabilidad, ganancia y pérdida de áreas (i.e., número de pixeles). Ya tenemos los mapas binarios para el periodo actual y el futuro (0, 1). Usaremos una función del paquete Biomod2 (BIOMOD_RangeSize).
bin.e1<-raster(bin.e1)
bin.e1.fut<-raster(bin.e1.fut)
bin.e1 <- reclassify(bin.e1, c(-Inf, 0, 0, 0, Inf, 1))
bin.e1.fut <- reclassify(bin.e1.fut, c(-Inf, 0, 0, 0, Inf, 1))
cambio <- BIOMOD_RangeSize(bin.e1, bin.e1.fut)##
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Range Size Computation -=-=-=-=-=-=-=-=-=-=-=-=-=
##
## Comparing 'proj.current' and 'proj.future'.
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=cambio## $Compt.By.Models
## Loss Stable0 Stable1 Gain PercLoss PercGain
## ensemble_weighted 30411 1301587 130205 162337 18.934 101.071
## SpeciesRangeChange CurrentRangeSize FutureRangeSize.NoDisp
## ensemble_weighted 82.138 160616 130205
## FutureRangeSize.FullDisp
## ensemble_weighted 292542
##
## $Diff.By.Pixel
## class : SpatRaster
## dimensions : 1363, 2741, 1 (nrow, ncol, nlyr)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -166.8333, -52.625, 14.54167, 71.33333 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : ensemble_weighted
## min value : -2
## max value : 1Vamos a mostrar el mapa.
Valores de -2 son los pixeles que se proyectan por perderse para la especie en el futuro. Valores de -1 son los pixeles que se proyectan por estar estables (no cambian). Valores de 0 son los pixeles que no ocupados actualmente pero tampoco en el futuro. Valores de 1 son los pixeles no ocupados actualmente pero proyectados por estarlo en el futuro.
plot(cambio$Diff.By.Pixel)
Hay varios paquetes para visualizar rasters.
https://oscarperpinan.github.io/rastervis/
levelplot(cambio$Diff.By.Pixel, par.settings=GrTheme())
TAREA:
Como tarea les queda hacer las proyecciones futuras con cada modelo invididual y evaluar qué tanta incertidumbre hay en cada algoritmo.
Eso es todo por ahora!